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Abstract 

Recent molecular dynamics simulation methods have enabled thermal conductivity of bulk ma- 
terials to be estimated. In these simulations, periodic boundary conditions are used to extend the 
system dimensions to the thermodynamic limit. Such a strategy cannot be used for nanostructures 
with finite dimensions which are typically much larger than it is possible to simulate directly. To 
bridge the length scales between the simulated and the actual nanostructures, we perform large- 
scale molecular dynamics calculations of thermal conductivities at different system dimensions to 
examine a recently developed conductivity vs. dimension scaling theory for both film and wire 
configurations. We demonstrate that by an appropriate application of the scaling law, reliable in- 
terpolations can be used to accurately predict thermal conductivity of films and wires as a function 
of film thickness or wire radius at realistic length scales from molecular dynamics simulations. We 
apply this method to predict thermal conductivities for GaN wurtzite nanostructures. 
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I. INTRODUCTION 



Thermal conducting properties of semiconductor nanostructures (e.g., nanowires) has 
been actively explored in recent years [IH5] because they directly impact many important 
applications including microelectronics and thermoelectrics. In the microelectronics appli- 
cation, a continuous decrease in feature sizes has resulted in a continuous increase in heat 
generation. This trend has placed an increasingly demanding requirement for the semicon- 
ductor materials to have a high thermal conductivity to effectively dissipate the excessive 
heat to the surrounding environment pQ . In the thermoelectrics [B] application, on the other 
hand, a low thermal conductivity is desired because it results in an increase in the energy 
conversion efficiency At the nanometer length scale, the effective thermal conductivity 
becomes sensitive to feature dimensions and defect concentrations. While this provides an 
effective means to tailor thermal conductivity for specific applications, the scaling of thermal 
conductivity against feature dimension is not always clear. Because experimental measure- 
ment of thermal conductivity is increasingly more challenging as the feature size decreases, 
a theoretical understanding of thermal conductivity as a function dimension can play a criti- 
cal role towards optimizing many nanostructure applications including microelectronics and 
thermoelectrics . 

To study the dimension effects on thermal conducting properties of nanostructures, pre- 
vious works have used the solution of Boltzmann partial differential equations [7HT2]. This 
approach is complex, requiring certain assumptions to reach simple analytical solutions. For 
example, the simple analytical equations for thermal conductivity provided in references 
[7H9] are only applicable either in a small or a large dimension, and hence they cannot be 
used to extrapolate data obtained from one dimension range to another. In addition, the 
Boltzmann partial differential equations involve certain input parameters, such as surface 
specularity, which may not be always available for a given material of interest. 

When an accurate interatomic potential is available, the use of molecular dynamics (MD) 
simulations in studying the thermal transport properties of crystals [T3H3T] may become de- 
sired. This is because the computational system used in MD simulations captures exactly 
the lattice nature of the crystal, which enables effects of surfaces and defects to be accu- 
rately incorporated. It has been shown that a reasonably accurate determination of thermal 
conductivity requires a real time of MD simulation for at least tens of nanoseconds |31j. At 
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this time scale, the system size that can be effectively employed usually contains no more 
than a million of atoms even with massively-parallelized MD simulations. For GaN, this 
translates to about 1.0 x 10 4 nm 3 of material volume. However, the GaN nanowires grown in 
experiments can have radius exceeding 20nm and length exceeding 20/im (32]. This corre- 
sponds approximately to a material volume exceeding 2.5 x 10 7 nm 3 . As a result, significant 
length-scale difference exists between the simulated and the real systems. 

Even with rather small systems, MD simulations have been relatively successfully applied 
to determine thermal conductivities of bulk materials based upon either the Green-Kubo 
(and its variations) [T3"HT§] or the "direct method" [20-31J. In the Green-Kubo method, 
periodic boundary conditions are used in all the three coordinate directions. As was demon- 
strated in previous work [31] and will be reexamined in the following, the use of periodic 
boundary conditions effectively extends the dimension of a small computational system to 
infinity. As a result, an infinitely large bulk crystal can be well captured with the Green- 
Kubo method even when a small simulated system is used. In the "direct method" , periodic 
boundary conditions can also be applied in the two coordinate directions perpendicular to 
the heat flux so that these two directions can be viewed as infinity. However, the "direct 
method" involves a heat source and heat sink along the heat conducting direction. A finite 
spacing, L, must be imposed between the source and the sink. Fortunately, both experi- 
ments and theories indicated that the inverse of thermal conductivity 1/k and the inverse 
of length scale 1/L satisfy accurately a linear scaling relationship [2"TW2"3"t |27J US] : 

i = i + f (l) 

K Kb L 

where k& is the bulk thermal conductivity at L — > oo , and a is a dimension-independent 
coefficient. To obtain bulk thermal conductivity, several simulations for different small cell 
lengths are performed, and the results can then be relatively accurately extrapolated to the 
infinite-size limit due to the linearity of the relationship. 

Periodic boundary conditions cannot be used for finite system dimensions. However, if a 
reliable linear scaling law that is applicable from nano up to macro scales is known, the ther- 
mal conductivity of finite systems at realistic length scales can still be accurately predicted 
based upon data obtained from MD simulations on a nano scale. An underlying assumption 
of Eq. Q is that the two dimensions perpendicular to the heat flux are infinite. Because 
of this, Eq. ([TJ is essentially a scaling law for 2D films where the heat flux is through the 



film thickness L. Unfortunately, no similar scaling laws were previously available for other 
heat flux directions (e.g., parallel to the film surface) or other nanostructures (e.g., nanowire 
or nanoparticles). As a result, previous MD simulations had not been applied to calculate 
thermal conductivities at realistic length scales for many interesting nanostructures, includ- 
ing cases where heat conduction occurs in the plane of a film or through the axis of a wire 

[Hoi Eg. 

Recently, we developed a theoretical scaling law that defines thermal conductivity of a 
nanostructure as a function of all of its three independent dimensions: thickness t, width 
W, and length L 
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where Kq jC , K\ fi , k 2jC , So, Si, S 2 , and d are seven constants that can in principle be determined 
from available thermal conductivity vs. dimension data. By performing very large MD 
simulations at different system dimensions, we have demonstrated that Eq. ^ is highly 
accurate from a nano scale all the way to macro scales [36J. The development of such an 
analytical scaling law has begun to enable MD simulations to be used to predict thermal 
conductivity of nanostructures at realistic length scales. The goal of the present work is 
threefold. First, we provide more detailed physics of the scaling law by adapting it for 
general 2D film and ID wire cases. Next, we explore the conditions and parameter space 
under which the scaling law can be accurately applied, and discuss the methods to predict 
thermal conductivity of nanostructures at realistic length scales. Finally, we perform large 
scale MD simulations to determine the [0001] thermal conductivities of a wurtzite GaN 
crystal constructed in two nanostructure configurations: (i) (1100) film with varying film 
thickness and (ii) a [0001] hexagonal wire with varying wire radius. GaN is chosen for the 
case study because it has excellent optoelectronic properties and can be easily integrated 
with the existing silicon structures. In addition, some GaN applications, such as laser diodes 
and high electron mobility transistors [37H42J. operate at high current and power densities. 
Understanding thermal transport of GaN nanostructures helps these applications. 
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II. SCALING LAW 



The underlying assumption of our scaling law |36j is accurate when the dimension of the 
structure is larger than the phonon mean free path. For GaN bulk crystals, the phonon 
mean free path has been estimated to be approximately 500 A at 300 K and 100 A at 500 K 
from both experimental data and kinetic theory jl3] . In nanostructures, the apparent mean 
free path is reduced due to surface scattering, causing the thermal conductivity to reduce. 
Interestingly, for given cross section, simple theoretical analysis indicated that the inverse 
of the phonon mean free path along the length direction is a linear function of inverse of 
the length [27]. This suggests that the inverse of the thermal conductivity is also a linear 
function of the inverse of the length, matching exactly the prediction of the scaling law. This 
means that the scaling law can actually be applied even when the length scale is comparable 
with the phonon mean free path. We will reexamine this in the following. 

Our theory can be extended and applied to arbitrary heat flux directions with respect 
to arbitrary shapes of the nanostructure. Here we confine our discussion to film and wire 
cases. 

A. The film case 

The geometry of the film case is illustrated in Fig. [IJa) where heat is assumed to flow 
through a finite length L of a box-shaped sample with a finite thickness t and an infinite 
width W — > oo. Note that for a true 2D film, L — > oo. A more general scenario of finite L 
is assumed here so that the theory can be applied with the direct method MD simulations 
where a finite spacing between heat source and heat sink must be used. It is recognized 
that the size-effect on thermal conductivity origins from the surface scattering of phonons. 
Hence, we separately consider surface and bulk regions of the sample. As shown in Fig. 
[jja), the sample is divided along the thickness direction into three smaller box-shaped 
regions (referred to as plates hereafter): the inner (core) plate has a thickness of t — 2d and 
is marked as "0" because it does not bound any y- surfaces, and the two outer (shell) plates 
have a thickness d and are designated as "1" because they bound one y- surface. When the 
thickness t is very large, we can always choose a sufficiently large shell thickness d so that 
the thermal transport behavior of plate is independent of the presence of the top and the 
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(a) core-shell plates of a sample 

low temperature T] 




W (-»«>) 



FIG. 1: Heat conduction through a finite length L of a box-shaped film with a finite thickness t 
and an infinite width W — > oo. 

bottom free surfaces that are far away. This means that the thermal conductivity of plate 
is independent of t and therefore can be expressed as a function of L only: k,q (L). The local 
thermal conductivity inside plate 1 is nonuniform near the surface. However, plate 1 still 
exhibits an apparent overall thermal conductivity. Note that at a large d, there is really no 
"distinguishable" interface between plates and 1 as the thermal transport properties from 
both sides of the interface approach the same bulk values. This means that once a large 
value of d is given, the apparent thermal conductivity of plate 1 can also be expressed as a 
function of L only: K\ (L). 

Using Fig. pT a), we assume that the left hand side of the sample is kept at a high 



6 



temperature of Th and the right hand side at a low temperature of TJ. Because the vertical 
temperature gradient is zero at the "indistinguishable" interface between plates and 1, we 
can list separately the thermal transport equations for the two types of plates: 

K o (L) = j^zf[ ( 3 ) 

L 

«x (L) = ^ (4) 

L 

where J and J\ are, respectively, the heat fluxes through plates and plate 1. Note that 
because of an assumed zero vertical temperature gradient at the 0/1 interface, the high and 
low temperatures are the same for both types of plates. The overall thermal conductivity 
of the system is expressed as 



J 



K 



Th—Ti 
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where the total flux J can be calculated as an area-weighted average 



W ■ it - 2d) ■ J + 2W ■ d ■ Ji _ it - 2d) ■ J + 2d ■ J x 
J ~ W^t " t (6) 

as in parallel conductors. Substituting Eqs. (|3]), Q, and ^ into Eq. (J5|, we have 

2d 

k (t, L) = n (L) - [kq (L) - K\ (L)) ■ — (7) 

Now we consider the thermal transport through the i-th plate (i = 0,1). Imagine that 
the plate is divided along the length direction into three sections: the center section contains 
a length of L — 2u and is marked as "c", and the two end sections contain a length of u 
and are marked shown in Fig. [l](b). Just as a subsurface thickness d subsumes 

the scattering of the side surfaces, a subsurface length u subsumes the scattering of the 
end surfaces (including the artificial effects of the thermostats). Here we distinguish d 
and u for generality. It can be seen that for a given large value of u, which is always 
possible when L is sufficiently large, the thermal transport behavior of the center section 
is independent of the presence of the two end surfaces that are far away. This means that 
the apparent thermal conductivity of the center section is equal to a constant /tj jC (i = 0, 1). 
In particular, kq >c corresponds to the bulk thermal conductivity Kb by definition. Similar 
to the discussion in the above, the apparent thermal conductivity exhibited by the two end 
sections is also independent of L and therefore is equal to another constant Because 
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heat flows through the three sections of the plate in serial, the heat flux J is a constant. We 
can list the temperature difference between the left and right ends as: 



Th - T, 



— • (L - 2u) + — ■ (2oj) 



(8) 



The inverse of the overall thermal conductivity of the plate, k { 1 (L), equals J 1 -(T/ l — T{) / L. 
We can therefore write: 

1 1 9m / 1 1 \ 
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Eq. ([£]) can be rewritten as 



where 
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(11) 



8{ combines the relative change of thermal conductivities between the center and the end 
sections with the length u, it therefore reduces one parameter. This reduction in parameters 



is expected because u and dependent. It can be seen from Eq. (11 lh that Si can be 



viewed as a characteristic length measuring the scattering of the end surfaces. Substituting 
Eq. (10) into Eq. (J7|), we have a scaling law for the thin film: 



K(t,L) 



L + 5 



2d 
t 



(12) 



L + S L + 5 X 

Eq. (12) is consistent with the previous work [51] . For instance, it reduces to Eq. Q 
when t — > 00, and it matches Eq. ^ when W — > 00. Eq. (12) involves five parameters 
ft^o.c «i,o <5cb <5i; and d. These five parameters have physical meanings and must be subject 
to some physical constraints. First, the surface thickness d and the associated surface 
thermal conductivity ki jC are dependent parameters. Hence, d can be selected and only the 
corresponding K\ yC value be treated as an unknown parameter. However, d is not completely 
arbitrary as it must be large enough to subsume the surface scattering effect. Once d is 



large enough, Eq. (12) can always predict accurate results regardless of its particular value 



for any film thickness t that satisfies the geometry condition t > 2d. On the other hand, a 



large d prevents Eq. (12) from being used for small thickness t due to the constraint t > 2d. 



So it is important to understand the low bound of d. Clearly d is sufficiently big if it equals 
the phonon mean free path in the bulk crystal. As described above, this might be overly 
stringent. 
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Once d is chosen, the remaining parameters can be fitted to the available data. For MD 
applications, it is necessary to perform several simulations at different dimensions in order 



to fit Eq. (12). Note that if the minimum system thickness used in these simulations is 
t m i n , then the largest d that still enables all the MD data to satisfy the geometry condition 
is t m i n /2. In order to find a small d to enable study of small structures, a trial-and-error 



approach can be used. For instance, d can be first set to t m i n /2, and Eq. (12) fitted to all 
MD data. If satisfactory fitting is obtained as will be described below, then the selected d 
is good. Otherwise d can be set according to the next thinnest sample, and the thinnest 
(tmin) sample is disqualified from the fitting. This process is continued until an appropriate 
d is found. There are also some useful relations. Because the end section is assumed to 
have more surface scattering than the center section, and plate 1 has surface scattering that 
is assumed to be insignificant in plate 0, we always have 5i > (i = 0, 1), K\ fi < «o,o 
and Ki (L) < k (L) (for any L). These conditions are automatically satisfied during fitting 
provided that the data to be fit are accurate and d satisfies the geometry constraint. 



Eq. (12) can be used for infinite 2D films. When t — > oo, the problem is essentially the 
heat conduction through the length L of a film (L is in fact the "thickness" in this case). 



When L — > oo, Eq. (12) gives thermal conductivity in the plane of a film as a function of 



the film thickness t. In addition, Eq. (12) can also be used for quasi- 2D cases or even ID 
cases, e.g. out-of-plane conduction, to explore the dimensional effects by using different t/L 
ratios. 



B. The wire case 

The wire case is illustrated in Fig. [2] where heat is assumed to flow through the length 
L of a circular sample with a finite radius r. Using the same theory described above, the 
sample is divided along radius direction into an inner, smaller cylindrical core with a radius 
of r — d and an outer cylindrical shell with a thickness of d. The core does not bound any free 
surfaces whereas the shell terminates with a radial surface, and hence the core and shell are 
denoted as "0" and "1" respectively. It can be seen that when r is very large, we can always 
choose a sufficiently large d to subsume the surface scattering effect so that the thermal 
conductivities of the core and the shell are independent of the wire radius and hence can be 
expressed as functions of L using kq (L) and K\ (L) respectively. In Fig. [2l we again assume 
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low temperature T[ 



high temperature T n 




FIG. 2: Heat conduction through a finite length L of a circular wire with a finite radius r. 

that the sample is held at a high temperature of at the left and a low temperature of T\ 
at the right. The thermal transport equations for the core, shell and overall system can be 
represented by Eqs. (j3]) - The total flux J, however, is modified as 

J • • (r - df + Ji • [tt • r 2 - tt • (r - d) 2 ] J • (r - d) 2 + Ji ■ d (2r - d) 



J 



7r • 



Substituting Eqs. (]3]), (]!]), and (13) into Eq. (15]), we have 

K (r, I) = n (L) - [«q (If) - «i (I)] • 



2d d 2 



(13) 



(14) 



where k (I) and «i (I) can be described by Eq. (10). Substituting Eq. (10) into Eq. (14) 
we have a scaling law for the wire: 

I ' K 0,c ( L ■ Kq c L ■ K\ 



K (r, I) 



(2d d 2 



(15) 



I + 5 \L + 5 L + 5 t J \r 
Eq. (15) is also consistent with the previous work [21] as it reduces to Eq. ([T]) when 
r — > oo, and it matches Eq. ^ using the geometry conditions of a circular wire [36J: 
2r, K2, c = ^i,c, and 82 = S\. Eq. (15) involves the same five parameters as 



t 



W 



in the film case. The geometry of the wire case, however, requires that d < r. For MD 
applications, the maximum d enabling all MD data equals the minimum radius r m j n used in 
the series of MD simulations. With d determined similarly as in the film case, the remaining 
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four parameters can be fitted to the measurements. If done correctly, the parameters should 
satisfy 5, > (i = 0, 1), /ti iC < k 0iC , and k x (L) < k q (L) (for any L). 



Eq. (15) can have numerous uses. When L — > oo, Eq. (15) represents thermal con- 



ductivity through an infinite ID wire as a function of wire radius. In particular, Eq. (15) 
indicates that thermal conductivity of wires is a linear function of 2d/r — d 2 /r 2 . When r is 
large, the thermal conductivity increases to a first order with ~ — 1/r, in agreement with the 



approximate equation derived by Lu et al from the Boltzmann equation [HE]- Eq. (15) can 
also be used in other cases. For instance, at r — > oo, the problem reduces to heat conduction 
through the thickness L of an infinite 2D film. It can be used for quasi- 2D films or even 
3D particles to explore the dimension effects by using different r/L ratios. 



III. MOLECULAR DYNAMICS METHODS 



One ultimate goal of our work is to enable MD simulations to predict thermal conduc- 
tivities of GaN films and wires at realistic, device length scales on the order of 100 — 1000 
A or more, which is not at bulk limit but too long to directly simulate with MD. Here we 
describe details of the interatomic potential used in the MD, the computational cells for film 
and wire configurations, and the thermal transport simulation method. 



A. Interatomic potential 

The previous work [SB EH] applied the Stillinger- Weber (SW) potential parameterized 
by Bere and Serra jSJ US] to calculate the thermal conductivity of GaN bulk crystals. To 
compare with the previous results, we also use the same potential in the present study. This 
potential gives reasonable prediction on dispersion relations, vibrational density of states 
(DOS), and heat capacity for bulk systems |31j . 



B. Computational system for films 

The computational system used for the film simulations is shown in Fig. |3j where the 
color scheme shows the temperature (red means the highest temperature and blue means the 
lowest temperature). Similar to Fig. [TJ we assume that the sample has a finite thickness t in 
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FIG. 3: Atomistic configuration for film MD simulations. 

the y- direction and an infinite width W — > oo in the z- direction. Following the customary 
approach with the direct method MD simulations [2"0TI3"T] . a periodic boundary condition 
was used in the x- direction. As will be shown below, this means that the system dimension 
in the x- direction is 2L rather than L. 

The equilibrium GaN has a wurtzite hexagonal crystal structure. The hexagonal crystal 
has three orthogonal directions [0001], [1100], and [1120]. To study thermal conduction 
along the [0001] direction of a low energy (1100) film, the computational supercell is aligned 
so that the x—, y—, and z— coordinates correspond respectively to [0001], [1 100], and [1120] 
directions. The experimental lattice constants of the hexagonal wurtzite are a = 3.19 A, 
c = 5.19 A, and internal displacement between Ga and N sublattices u = 0.377 [IB]. With 
the SW interatomic potential used here, the zero temperature lattice constants are a = 
3.19 A, c = 5.20 A, and u = 0.375. Converting the hexagonal crystal to the smallest 
orthogonal unit cell, the lattice constants of the unit cell are respectively a\ — c — 5.2000 A, 
a2 = 2 • a ■ cos (vr/6) = 5.5252 A, and = a = 3.1900 A in the x— , y—, and z— directions. 
For convenience, the system dimension can be represented by the number of cells n\, n,2, 
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and 713 in the x—, y—, and z— directions. In addition to the difference in unit, rii, 712, and 
n 3 always refer to the simulated size whereas t, W and L refer to the real size that can 
become infinite. Two series of sample dimensions, one corresponding to n\ = 100 and the 
other one corresponding to rii = 150, were studied at various i%2 ranging from 50 to 150 
and in addition at t — > 00 using a fixed = 5. The film scenario with the (1100) film 
surfaces was simulated using a free boundary condition in the y— direction with surfaces 
terminated between the larger spacing as will be described in details in the wire case. Such 
a termination ensures stable surfaces and therefore no surface reconstruction was observed 
in simulations as supported by the experiment [47J. The t — > 00 or W — > 00 case was 
simulated by using a periodic boundary condition in the corresponding direction. Note that 
although the periodic boundary condition is also used in the x— direction, the meaningful 
dimension in the x— direction for the direct method MD simulations is the spacing between 
heat source and heat sink. This spacing is not extended by the periodic boundary condition 
and is always finite. 

C. Computational system for wires 

Unlike the circular wire assumed in Fig. [2j the GaN wires observed in experiments are 
often hexagonal with (0001) axis and {1100} facets. The computational system used for 
simulations of such wires is shown in Fig. (4^a). The crystal orientations are the same as 
those used for the film simulations, except that atoms beyond specified {1100} surfaces are 
removed. As in Fig. [2j the system is assumed to have a finite radius r, which is defined 
as the minimum distance between the center of the wire and the surface, Fig. |4](b) . Again 
the system dimension in the x- direction is assumed to be 2L rather than L to facilitate the 
periodic boundary condition. 

The cross section of the hexagonal wire is examined in details in Fig. |4^b) . It can be seen 
that the {1100} atomic planes have two different spacings. We found that if the surface is 
terminated between the small spacing, such as that shown by the dash line, then the surface 
is composed of many dangling bonds resulting in high energy and unstable configurations. 
Hence, the surfaces of our wires are always terminated between the large spacing. It can 
also be seen that the smallest repeatable spacing of the {1100} planes is a r = a ■ cos (tt/6) = 
2.7616 A. For convenience, the system dimension in the x— and the radial directions can be 
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(a) nanowire geometry 




^> z[1120] 
(b) cross section 

/\\ y [Tioo] 



unstable surface if cut here 



small 




FIG. 4: Atomistic configuration for wire MD simulations (black and white colors distinguish Ga 
and N atoms, and the star shows the center of the wire cross section). 

respectively represented by number of cells n\ and n r (in unit of a r ). A matrix of dimensions 
with the longitudinal dimension n\ ranging from 100 to 300 and the radial dimension n r 
ranging from 4 to 12 was explored. Here free boundary conditions were used in the y— and 
z— directions and the periodic boundary condition used in the x— direction. 



D. Heat transport simulation algorithm 

The thermal transport MD simulations were performed under a constant number of 
atoms, constant system volume, and constant system energy (NVE) condition using a time 
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step size of At = lfs. To accurately account for the effect of thermal expansion and elim- 
inate the errors due to statistical fluctuation of the simulated temperature, the following 
steps were used to create the initial crystal. First, a crystal was created by assigning atom 
positions according to the prescribed crystal lattice and the known lattice constants at 
zero temperature. A molecular dynamics simulation in the constant atom number, (zero) 
pressure, and temperature (NPT) ensemble was subsequently performed for a total of 20ps 
period. The desired simulated temperature was achieved using the velocity rescaling method. 
After discarding the first lOps simulation to allow the system to reach a steady state, the 
average crystal sizes and average total (kinetic and potential) system energy were then cal- 
culated for the remaining lOps. We then created another crystal using the average sizes 
obtained at the finite temperature. We can also calculate the potential energy of this newly 
created crystal by simply performing an energy calculation simulation. The difference be- 
tween the average total system energy obtained from the NPT run and the potential energy 
of this crystal prescribes exactly the amount of the kinetic energy that needs to be added 
in order for this crystal to exhibit the same total average energy over the subsequent long 
constant energy thermal transport simulation. We added precisely this amount of kinetic 
energy into the system by first assigning velocities to atoms according to Boltzmann dis- 
tribution and then rescaling the velocities under the zero total linear momentum condition 
[26l HU H9] . The thermal transport simulation is started immediately without the conven- 
tional long NPT or NVT simulation to establish the initial temperature. The advantage of 
this approach is that once steady-state is reached, the average temperature of the system 
matches exactly the desired temperature. In practice, we found that the difference between 
the average temperature in the middle of the heat source and heat sink (where the thermal 
conductivity was calculated) and the desired temperature is well below 1 K (often near 0.01 
K or less). This method has the same accuracy level as the "doubling temperature" method 
applied previously [31J. It is more general because the initial crystal is not required to be 
in the minimum potential energy configuration, which may not be easily determined when 
the system includes surfaces. 

The direct method requires the creation of a heat source and a heat sink. As shown in 
Fig. |3j the heat source corresponds to the red region at the far left of the cell, and the heat 
sink is the blue region near the middle of the cell. During simulations, the size and location 
of the source and sink regions are specified. With appropriate choices of system length, 
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source and sink region width, and locations of different regions, we ensured that the source 
and sink regions were geometrically identical, and that the left side of the sink (or source) 
region was exactly symmetric to its right side, up to another source (or sink) region (which 
may be its periodic image). Here the width of the source or sink region is around 40 A. It 
can be seen from Fig. [3] that the heat source at the left side has an image at the right side 
under the periodic boundary condition in the x- direction. Consequently, even the length of 
our system is 2L, the spacing between the source and the sink is still L as shown in Fig. [TJ 

The constant flux [2"5T - |2l?| |4*S] method was used to create the temperature gradient. In this 
method, a constant amount of energy is added to the hot region and exactly the same amount 
of energy is removed from the cold region at each MD time step using velocity rescaling 
(while preserving linear momentum). To ensure that the high and the low temperatures are 
reasonable and consistent in different runs, the heat flux has been adjusted within 0.00035 
to 0.0010 eV/(ps ■ A 2 ) to give a consistent high-low temperature difference (say 5-10 K). 
To generate extremely accurate results, the first 0.4 ns simulation was discarded to allow 
the system to reach a steady state, and the remaining duration of simulations was chosen 
to be at least 11 ns and many reached over 20 ns. To compute the temperature profile, the 
system dimension in the x— direction is divided into a grid. The temperature of each of 
the grid cells was averaged over the remaining time of simulations. The temperature profile 
and the input heat flux were used to calculate the thermal conductivity using Fourier's Law, 
Eq. To estimate the statistical error of the calculated thermal conductivity, the total 
averaging time was divided into 20 subsections and thermal resistivity (or conductivity) was 
calculated for each of the subsections. From these data, analysis of statistical error was 
performed. For more details on the procedure please see [31 J. 

It should be noted that under the periodic boundary condition in the x- direction, the 
observed dependence of thermal conductivity upon the system length L comes primarily from 
the scattering of the interfaces at the hot and the cold regions. The functional dependence 



on L can still be well described by Eqs. (12) and (15), albeit 6i should be viewed as 
an interface scattering parameter rather than surface scattering parameter. Because our 
analysis extrapolates the MD data to an L — > oo limit (i.e., true film and true wire), the 
interface approximation will not affect the results. 
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TABLE I: Parameters describing the general scaling law. 



structure 


T(K) 


d(A) 


K , C {W/K ■ m) 


Ki, c (W/K ■ m) 


So (A) 


*i (A) 


film 


300 


138.13 


178.38 


151.65 


1288.10 


1329.65 


film 


500 


138.13 


75.42 


56.16 


623.38 


471.74 


wire 


500 


11.05 


75.42 


47.30 


615.51 


1324.45 


wire 


500 


22.09 


75.42 


50.30 


604.71 


795.87 



IV. RESULTS 



A. Film thermal conductivity as a function of film thickness 

Systematic MD simulations were performed to derive thermal conductivity of film as a 
function of film thickness at two temperatures of T = 300 K and T = 500 K. As the Debye 
temperature was estimated to be in the range 350-600 K [13J |50j EE] , the system is expected 
to behave classically especially at 500 K. Although 300 K is at the lower bound of the 
estimated Debye temperature range, it is chosen for study because the low temperature data 
has less thermal fluctuation and therefore can provide stronger model validation. At each 
temperature, two series of MD simulations corresponding respectively to 2L m 520A (n\ = 
100) and 2L rj 780A (n\ = 150), were performed at different thickness t but fixed width 
W — > oo as described in the above. Previous MD simulations have determined the thermal 
conductivities at various length L but a fixed thickness t — > oo and a fixed width W — )■ oo 



Both present and previous data are collectively used to fit Eq. (12) using a chosen value 
of d, and the results of the fitted parameters are shown in Table [Tj Both the MD data and the 
fitted curves are shown in Figs. [5] and [6] for the 300 K and 500 K temperatures respectively. 



Eq. (12) indicates that at a fixed L, thermal conductivity k is a linear function of 1/t, and 
at a fixed t, the inverse of thermal conductivity 1/k is approximately a linear function of 
1/L. Hence, Figs. [5][a) and |6ja) show the k vs. 1/t plots at fixed L whereas Figs. [5| 



and[6tb) show the 1/k vs. 1/L plots at fixed t. It can be seen that the linear relationships 



predicted by Eq. (12) is strikingly reproduced by the MD data, and the agreement between 
the fitted lines and the MD data is excellent using a single set of parameters (d, «o,o ^i,o 
8q, and 5i) for both the thickness and the length functions. 

In Fig. [5] and [6| the shaded area indicate the length-scale no-man's land, meaning that 
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FIG. 5: Thermal conductivity data obtained for film at 300 K. 
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FIG. 6: Thermal conductivity data obtained for film at 500 K. 

the system dimension is too large to be directly simulated using the MD methods. Fig. 
[5jb) and [6]^b) show that the linear scaling law can be used to predict thermal conductivity 
in the shaded region through extrapolation based upon the MD data obtained in a small 
dimension range. Although the linear extrapolation is expected to produce reliable results, 
the accuracy of the extrapolated values is difficult to confirm directly. In sharp contrast, 
Fig. [5]^a) and [6](a) show that because the thermal conductivity at an infinite thickness t 
can be obtained in the MD simulation by using periodic boundary condition and because 
the linear relationship extends to this value at t — > oo, an extremely reliable interpolation 
can be used to predict thermal conductivity at any thickness dimension. Now we can see 
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(a) 300 K (b) 500 K 




film thickness t (A) film thickness t (A) 



FIG. 7: Predicted thermal conductivity of film as a function of film thickness. 

that the satisfaction of the linear relationships by the MD data and the extrapolation of the 
linear relation to the bulk limit can be used to determine if the selected d value is sufficiently 
large. 

For true films, we set L — > oo. Using the parameters listed in Table [T] and Eq. (12), 



thermal conductivity of film was calculated as a function of film thickness, and the results 
obtained at 300 K and 500 K temperatures are shown respectively in Figs. [7[a) and[7^b). 

To confirm our hypothesis that the boundary scattering region is confined to the surface, 
we calculated the heat flux of a thin film in slabs parallel the surface of the film. The calcu- 
lated relative flux and temperature profiles are shown in Fig. [8] Clearly it can be seen that 
flux remains nearly constant in the interior but is degraded near surfaces. Temperature, on 
the other hand, is constant across the entire thickness of the sample. This strongly validates 
that the conductivity in the interior is a constant whereas it is reduced at surfaces. This 
phenomenon is unknown in the past. It indicates that although the thermal conductivity at 
small scale can be thought to depend on the phonon mean free path, it can still be mani- 
fested through the core-shell phenomenon as assumed in our model. This accounts for why 
our scaling law agrees well with the effects deduced from considering only the phonon mean 
free path [27] . 
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FIG. 8: Flux and temperature profiles along the thickness of a thin film. 



B. Wire thermal conductivity as a function of wire radius 



Having verified the film case at both 300 K and 500 K, MD simulations were performed 
to study thermal transport for wires at 500 K. Only one temperature is studied due to the 
expense of the calculations. 500 K was chosen because it is likely to be above the Debye 
temperature. To deduce all the parameters, simulations were carried out using a matrix of 
three sample length of 2L « 520A (m = 100), 1040A (m = 200), and 1560A (m = 300) and 
four sample radius of r ~ 11 A (n r = 4), 16A (n r = 6), 22A (n r = 8), and 33A (n r = 12). Eq. 



(15) indicates that at a given length L, k is a linear function of 2d/r — (d/r) and at a given 



radius r , 1/ k is approximately a linear function of 1/L. Here we explore two different chosen 
d values, d = 11.05A and d = 22.09A. Using d = 11.05A (which is the largest d that enables 
all MD data to satisfy the geometry condition of the scaling model), k vs 2d/r — (d/r) 2 
curves at different L are shown in Fig. [9]^a) and 1/k vs 1/L curves at different r are shown 
in Fig. |9^b). In addition, the k values obtained at infinite cross section dimension (i.e., 
r —7- oo) obtained in the previous work [31] are included in Fig. [9j In Fig. [9j the lines 
are calculated using Eq. (15) with fitted parameters displayed in Table [l| It can be seen 



that Fig. 9 exhibits some linear relationships predicted by Eq. (15). However, the overall 



match between the MD data and the model prediction is not great. Most seriously, Fig. 
9^a) indicates that a linear regression using merely the data points at large 2d/r — (d/r) 2 
values would not closely extrapolate to the data point at 2d/r — (d/r) 2 = (i.e., r — > oo), 



20 



(tt) effect of thickness 



55 

? 50 
X 45 

^-40 
X 

^35 
.jj 30 

c 

g 20 

1 15 
I 10 

5 
(I 





1 1 
shaded region: 


-5- L-520A 
o-L~1040A 




radius too big to simulate 

E 






-*-L~ 1560A 




L,. ■ 








j — , 












--•Oi, 




T = 500 K 











effect of length 

0,14 



S 0,12 

e, o.io 

0.08 

.« 0.06 
cn 

a 

"5 0.04 
■S 0.02 
0.00 



1 


^- length too big to simulate 


-o- r~ 11 A 
- -tt - r~ 16 A 
r~22A 
— «- r~28 A 
-o- r~33A 
— •— r co 


..-'**" ^^--3h=_ 




T = 500 K 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
unitless radius scale 2d/r - (d/r) 2 



0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040 
inverse of length 1/L (1/A) 



FIG. 9: Thermal conductivity data obtained for wire at 500 K, d = 11.05A 



and Fig. |9jb) shows significant deviation of the predicted curves from the data points. 

The problems exhibited by Fig. [9] is inherently related to an important condition: the 
parameter d must be chosen to be sufficiently large to subsume the surface scattering effect 



in order for Eq. (15) to be valid. To explore this further, we increased d to 22.09 A. This 
disqualified the data obtained at small radii of r fa 11 A and r fa 16 A. The remaining data 
were fitted to Eq. (15), and the parameters thus obtained are included in Table [T] With the 
increased value of d, the results similar to those in Figs. |9](a) and |9](b) are recalculated and 
are presented in Figs. [lO^a) and [To](b) respectively (note that the scales in the horizontal 



axis are different between Figs. [9] and 10 due to the use of different d values). Clearly, a 
significant improvement is achieved. In particular, Fig. [To^a) indicates that the data points 
obtained at large 2d/r — (d/r) 2 values can be linearly connected to accurately extrapolate 
to the point at 2d/r — (d/r) 2 = 0, and Fig. [To|b) shows an improved agreement between 
data points and the predicted curves. 

The value of d used in the wire case is constrained by the radius r used in the MD 
simulations. Unlike the film case where the thickness t is independent of the width W, 
increasing the radius in the wire case would increase the system dimensions in both y- and 
z- directions. As a result, the magnitude of the radius is more severely constrained by the 
computational expense. In general, current computer resources may not permit the use of 
radius significantly above n r = 12. The d value permitted by this radius is still small and 



significant further improvement over the one seen in Fig. 10 is likely to be achieved if d 
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FIG. 10: Thermal conductivity data obtained for wire at 500 K, d = 22.09A 

can be substantially larger. To improve the results, we explored an alternative method by 
again using the scaling theory. Film simulations described in the previous section resulted 
in the determination of five parameters Ko, c , ki iC , 5q, Si, and d. These five parameters have 
physical meanings and are invariant in the wire configurations. As a result, we can directly 



use these parameters and Eq. (15) to calculate wire thermal conductivity. Note that due 
to the relatively low computation cost, the d value used in the film case reaches 138 A, 
Table [T| As a result, significantly more accurate results can be expected. Some results of 

^b) using 
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this calculations at a few different radii r > d = 138.13 A are included in Fig. 
stars. It can be seen that although the predicted results from film simulations using big 
d = 138 A and wire simulations using small d = 22.09A are different, the difference is 
relatively small at least for GaN. Note that while the film parameters based upon a big 
d = 138 A give better results, they cannot be used to calculate thermal conductivity at 
r < d = 138 A. 



For true wires, we set L — > oo. The film parameters and Eq. (15) were used to predict 



wire thermal conductivity as a function of radius. The results are shown in Fig. 11 
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FIG. 11: Predicted thermal conductivity of wire as a function of wire radius. 
V. DISCUSSION 

A. Generalized scaling law 

The scaling law described in the above can be generalized to any prismatic wire with ar- 
bitrary cross-section. Approximated as coarse-grained resistance network with temperature 
being a function of x-coordinate only, the cross section of such a wire can be broken up into 
a finite number of distinct areas Aj, with one area Aq completely in the interior of the wire 
and all the others i > touching the boundary of the wire. The total cross-sectional area is 
then A = Ai. The total flux J is partitioned such that 



where the second equality is obtained from Fourier's law, and «y is the conductivity of the 
z-th area and the j-th section. Here we use j = -1 and 1 to indicate the two end sections 
and j = to indicate the middle section. It is assumed that the two end sections are similar 
so that Ki-i = k^i = Kj,±i and AT_i = ATi = AT±i. Under the geometry conditions that 
Ax_i = Ax\ = d and Ax-± + Axq + Ax\ = L, we have 




(16) 




(17) 
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Using the continuous temperature field condition AT = 2AT±i + AT , Eq. (17) can be 
written as 



2d 



L — 2d 



The constituent areas can be classified into three types: corners, boundaries, and the interior. 
Each relates to the overall radius R of the wire. In a self-similar fashion, the interior scales 
with R 2 ; the boundaries scale with R; and the corners are essentially constant. Using 
geometry-specific constants Cj, the areas of different types can be expressed as A^/A = Ci/R 2 
for the corners, Ai/A = Ci/R for the boundaries, and A /A = 1 — Y2i>o Ai/A for the interior. 
The overall thermal conductivity of the wire is then 
JL 



AT 



2d/L 



« ,±i(l - l/RT, ieB Ci - l/R 2 £ ieV G) + V^EieB «i,±iC< + l / R2 E ie v «i,±iC< 



+ 



2d/L 



(19) 



Kofl(l ~ 1/REieB Q - l/R 2 E, e v C t ) + 1/REteB KiflCi + l/R 2 Eiev ^oQ. 
where B and V represent the sets of boundaries and corners respectively, L and R are input 
sample dimensions, and the rest of the parameters, i.e. d, Ci, and Kij, need to be fitted. 



It can be seen that Eq. (19) correctly reduces to the bulk thermal conductivity value kqq 



as L — > oo and R — > oo. Practical nanowires usually have symmetric cross section, thus it 
is possible to reduce the number of independent parameters. For an equilateral triangular 
prismatic wire, for example, the free parameters would include Ko,±i an d ^o,o fo r the core, 
fci,±i and Ki t Q for the boundaries, and K2,±i and ^2,0 for the vertices in addition to d, C\, 



and C2. Eq. (19) then becomes 



K 



2d/L 



ko,±i(1 - C x /R - C 2 /R 2 ) + d ■ K Xt±1 /R + C 2 ■ « 2> ±i/i? 2 



+ 



2rf/L 



ko,o(1 - Ci/i2 - C 2 /i? 2 ) + d • K lt0 /R + C 2 ■ K 2fl /R 2 



(20) 



B. EfTect of periodic boundary conditions 

The use of periodic boundary condition eliminates surfaces. Strictly speaking, however, 
there is still an interface between the periodic image that is not exactly the same as if 
there were no such an image, e.g. phonon wavelength is restricted by the size of the periodic 
cell. The effect of this interface has been explored previously by using the periodic boundary 
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FIG. 12: Effect of cross section area on thermal conductivity under the periodic boundary condition. 



conditions in the two cross section directions with different simulated cross section areas [21] ■ 
The results showed that the periodic interfaces parallel to the heat flow have no significant 
effect on thermal conductivity [31 J . As a result, periodic boundary condition was thought to 
be able to extend the system dimension to infinite. The past study, however, only examined 
the simulated cross section area between 180 A 2 and 890 A 2 . In the present work, the cross 
section area reached as high as about 13250 A 2 . Hence, we re-examine the effect of cross 
section area under periodic boundary conditions. Thermal conductivities were calculated at 
300 K using periodic boundary conditions in both y- and z- directions with fixed simulated 
dimension of ri\ = 150 (2L ~ 780A), n% = 5 (~ 16 A) and various thickness of ri2 between 3 



(f=s 16A) and 150 (~ 830A). The results are shown in Fig. 12. Fig. 12 confirms that when 
the periodic boundary condition is used to extend the cross section dimension to infinity, 
the magnitude of the simulated dimension does not affect the thermal conductivity in a wide 
range of cross section between 260A 2 and 13220 A 2 . Furthermore, since we only varied the 
thickness t, our result provides a strong evidence that the aspect ratio t/W also does not 
affect the thermal conductivity under the periodic boundary condition. 
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C. Effect of surface stress 



The simulations described above are based upon a Stillinger- Weber potential. An essential 
feature of SW potential is that its pairwise functions decay to negligible values within the 
second nearest neighbor distance of atoms. A problem with the nearest-neighbor potentials 
is that they predict zero surface stresses [52] • The neglect of surface stresses may alter 
the thermal conductivities of nanostructures. We have also explored the use of alternative 
potentials such as the Tersoff GaN potential [53] in our studies. Unfortunately, we discovered 
that the Tersoff potential severely underestimates the thermal conductivity as compared 
with the experimental measurement. With the SW potential clearly the better choice for 
the thermal transport simulations, we carefully assess the surface stress effect. This was 
done by modifying the SW potential to artificially create a surface stress. All the pairwise 
functions used in the SW potential are expressed in terms of the scaled atomic spacing as 
F(R/a) [HI H5], where R is the atomic spacing and a is a scaling factor. As a result, we 
scaled the parameter a by 0.98 to reduce the equilibrium bond length by 2 %. We then used 
the modified potential for the two atomic layers of atoms on the nanostructure surfaces and 
the unmodified potential for the remaining atoms. This effectively created a surface tensile 
stress. This is a good approximation because it realistically captures the surface stress in 
the subsurface region despite the variation of atomic interaction of the surface atoms. This 
modified scheme of interatomic potential was then used to perform a selected case of film 
simulation at a temperature of 500 K, a sample length rii = 150 (2L w 780A), a width n 3 = 5 
(16 A), and a thickness ri2 = 100 (~ 550A). To understand the surface stress of the sample, 
the system was first relaxed using molecular statics energy minimization simulation at a 
constant volume determined from an NPT MD simulation at 300 K. The hydrostatic stress 
was calculated using the Virial theorem [51] and was binned along the thickness direction. 
The surface stress was calculated as the difference in hydrostatic stress between surface and 
bulk. This calculated surface stress is plotted as a function of position along the thickness in 
Fig. [13] for both modified and unmodified SW potentials. It can be seen that the modified 
potential clearly shows a significant surface stress compared with the unmodified potential. 
Thermal conductivities were calculated and we found 27.56 ± 0.36W/K ■ m for the modified 
potential as compared with 27.70 ± 0. 43Wy.fr ■ m for the unmodified potential. Based upon 
the data, we do not expect that neglecting the surface stress by the SW potential significantly 
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FIG. 13: Stress profile along the film thickness. 



affects the thermal conductivity estimates. 



VI. CONCLUSIONS 



We have explored general scaling equations that explicitly express thermal conductivity 
of film and wire as functions of dimensions. Based upon these scaling equations, we have 
demonstrated methods that enable molecular dynamics simulations to be used to predict 
thermal conductivity of nanostructures at realistic length scales. We have performed ex- 
tensive MD calculations of thermal conducting properties along the [0001] direction of GaN 
films and wires. The following conclusions have been obtained: 

• The linear relationships predicted from the scaling equations hold extremely well for 
MD data at a large parameter d. Reliable prediction of film thermal conductivity as 
a function of film thickness has been achieved using linear interpolation. 

• Thermal flux in nanostructures exhibit a clear difference between the surface and core 
whereas temperature is nearly constant across the cross-section, thereby verifying the 
core-shell assumption of the scaling theory and the near-zero heat flux perpendicular 
to the axis connecting hot and cold reservoirs. 
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• Due to the limitation of computational cost, parameters deduced from direct MD sim- 
ulations of wires may not sufficiently accurately predict the wire thermal conductivity 
at large wire radii. However, the parameters deduced from film simulations enable the 
derivation of a reliable expression of wire thermal conductivity as a function of wire 
radius. 

• The simulated dimension does not affect the thermal conductivity when the dimen- 
sion is transverse to the heat flow and a periodic boundary condition is used in that 
direction. Hence, the periodic boundary conditions can be used to accurately extend 
the system dimension to infinity. 

• The surface stress due to the contraction of surface bonds does not sensitively affect the 
thermal conductivity. As a result, the SW potential should be a sufficiently accurate 
force field for surface thermal transfer problems. 
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